Capillary waves at the liquid-vapor interface and the surface tension of water models 
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Capillary waves occurring at the liquid- vapor interface of water are studied using molecular dy- 
namics simulations. In addition, the surface tension, determined thermodynamically from the differ- 
ence in the normal and tangential pressure at the liquid-vapor interface, is compared for a number 
of standard three- and four-point water models. We study four three-point models (SPC/E, TIP3P, 
TIP3P-CHARMM, and TIP3P-Ew) and two four-point models (TIP4P and TIP4P-Ew). All of the 
models examined underestimate the surface tension; the TIP4P-Ew model comes closest to repro- 
ducing the experimental data. The surface tension can also be determined from the amplitude of 
capillary waves at the liquid-vapor interface by varying the surface area of the interface. The sur- 
face tensions determined from the amplitude of the logarithmic divergence of the capillary interfacial 
width and from the traditional thermodynamic method agree only if the density profile is fitted to 
an error function instead of a hyperbolic tangent function. 



I. INTRODUCTION 

The ability to derive accurate property predictions for 
the liquid-vapor interface is a key test for an atomistic 
force field. Because of the frequent occurrence of water 
in systems of chemical and biological interest, interfacial 
property prediction is especially vital for force fields of 
water. The most important of these properties is surface 
tension, an intensive quantity that measures the differen- 
tial surface work required to increase the interfacial area. 
Accurate models of the surface tension of water are essen- 
tial for conducting large-scale simulations of the wetting 
and spreading of water droplets at surfaces. 

An interface between two distinct thermodynamic 
phases can be characterized by a local gradient of an 
order parameter whose mean value changes between 
phases, such as the boundary between a liquid and its 
own vapor below the critical temperature T c . For sim- 
ple fluids, thermodynamic arguments predict that the in- 
terfacial width A depends only on temperature and the 
interaction energies within each phase and across the in- 
terface. However, the presence of the interface breaks the 
translational invariance of the system, inducing Gold- 
stone fluctuations or "capillary waves" at the interface 
P, • Density functional studies suggest that the sur- 
face tension predicted by capillary-wave theory exhibits 
a minimum as a function of the wavelength Q . Previous 
studies of capillary waves involving water have tended to 
focus on liquid- liquid interfaces or on model fluids 0, Q , 
and have generally examined relatively small systems of 
less than 10, 000 molecules; the present study represents 
the first study of capillary-wave behavior at the liquid- 
vapor interface of water. 

For two-dimensional interfaces, these noncritical fluc- 
tuations give rise to a logarithmic increase in the inter- 
facial width A with increasing L», the length of the in- 
terface. Most previous simulations 0, of the liquid- 
vapor interface in three dimensions did not investigate 
the dependence of A on the size of the interface. The 
purpose of this paper is to present atomistic molecular 
dynamics (MD) simulations of the liquid-vapor interface 



of water. In particular, we obtain the surface tension 7 
in two different ways: from the difference in pressure par- 
allel pu and perpendicular p± to the interface {"f p ), and 
from from the dependence of A on Lu (7™). We confirm 
the previous result that j w depends on the functional 
form chosen to fit the order parameter (density profile) 
through the interface @ . In particular, fitting the order 
parameter to an error function gives results for j w which 
are in strong agreement with j p . However, fitting our 
data to a hyperbolic tangent function, a functional form 
derived from mean- field arguments 0, gives results for 
7u, which are systematically smaller than j p and further 
away from experimental results. 

There are currently a large number of different atomic 
models for water. Guillot provides an extensive list of 
models developed through 2001 0; several additional 
models have been introduced since then 0, 0, 0] . 
The simplest of the commonly-used atomic models, the 
SPC model is a rigid three-point model with fixed 
charges; the most complex model, the POL5 model fl5| . 
is a polarizable five-point model. Most of the commonly 
used models are three- and four-point models. In three- 
point models, such as SPC/E and TIP3P [13 , the 
electric charges are assigned directly to the hydrogen and 
oxygen atoms; four-point models, such as the TIP4P 
|17| and Watanabe-Klein [18| models, locate the negative 
charge at a massless point a fixed distance away from the 
oxygen atom. Five-point models, such as TIP5P [l9| . 
and the early Bernal-Fowler [2(j and ST2 [2l| models, 
represent the negative charge of the oxygen using a pair 
of massless charges to capture the quadrupolar behavior 
of water. Polarizable models, including the SPC/FQ and 
TIP4P /FQ models 22] , allow the magnitude of the point 
charges to be treated as variables which can fluctuate ac- 
cording to the local environment. 

The proliferation of models has been motivated largely 
by the need to reproduce various physical and thermody- 
namic properties, such as the bulk density, the oxygen- 
oxygen radial distribution function, the heat of vaporiza- 
tion, and the diffusion coefficient. However, some models, 
such as the recent TIP3P-Ew and TIP4P-Ew models 
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[l2| , are reparameterizations of existing models designed 
to account for changes in the treatment of long-ranged 
electrostatic interactions. 

Most of the available water models adequately repre- 
sent at least some of the thermodynamic properties of 
water; for a comprehensive review, see Jorgensen et al. 
[l3|. Kuo et al. have shown that the changes intro- 
duced between, for example, the TIP4P and TIP4P/FQ 
models have little influence on properties such as the 
bulk liquid density or the mean distance between oxygen 
atoms either in bulk or at the interface ^23]. However, 
the simulation behavior of models with nearly identical 
parameters can be markedly different: Mark and Nils- 
son have noted significant variation in physical and ther- 
modynamic properties such as the self-diffusion constant 
and the radial distribution function of various three-point 
water models [3 H3 • Less is known about how well the 
various water models describe the liquid-vapor interface 
and the surface tension 7. Our own work, however, sug- 
gests that even models with very similar density and dis- 
tribution profiles can have quite different predictions for 
surface tension. 

Experimental studies demonstrate that the surface 
tension of water decreases with a slight quadratic depen- 
dence on temperature in the range 273 K < T < 373 K 
IH EllH lH H3. Surface tension results for higher 
temperatures have not been reported in the literature; 
we extrapolate the reported experimental data to higher 
temperatures. There have been a few studies for vari- 
ous three-point models [l^mmmmillllilll 
which show that while most water models reproduce the 
observed decrease in surface tension as temperature in- 
creases, they tend to underestimate 7 by amounts be- 
tween 25 and 50 percent. Only Alejandre et at, Huang et 
al, and Shi et al. £nim>IH3 report adequate agreement 
with experimental data. However, as we show below, the 
apparent agreement of both Alejandre et al. [3lJ and 
Shi et al. [37J is the result of inadequate simulation time. 
Alejandre et al. also employ a reciprocal-space mesh that 
is too coarse, while Huang et al. report values only for 
the SPC and SPC/E models at 298 K H3. 

Our primary goal is to study capillary waves at the 
liquid-vapor interface of water, and to distinguish be- 
tween various functional representations for the density 
profile near the interface. Additionally, we first deter- 
mine the surface tension as a function of temperature for 
six commonly used three- and four-point models of wa- 
ter, in part to establish a basis for comparison with the 
capillary- wave simulations. 

In Section II, we provide a brief overview of methods 
for computing the surface tension from molecular simu- 
lation data, of the various water models examined in this 
study, and of the simulation methods employed. Section 
III presents our findings on the temperature dependence 
of the surface tension, as well as the effects of the tail 
correction, interaction cutoffs, and reciprocal-space mesh 
refinement. We discuss the results obtained from the 
analysis of capillary waves at the liquid-vapor interface 




FIG. 1: (Color online) Sample simulation cell, showing equi- 
librated configurations of 1000 SPC/E water molecules at 
(a) 300 K and (b) 500 K. The dimensions of the cell are 
L x = L v = L11 = 2.3 nm and L z = L± — 13.5 nm. 



in Section IV before offering our conclusions in Section 
V. 



II. MODELS AND METHODOLOGY 
A. Surface tension 

1. Thermodynamic method 

There are two primary methods used to compute the 
surface tension using molecular simulations. The first ap- 
proach, developed by Tolman [3fll | and refined by Kirk- 
wood and Buff [40j. computes the surface tension as an 
integral of the difference between the normal and tangen- 
tial pressures p± (z) and p\\ (z): 



1p 



(p± ( z ) ~P\\ 0)) dz, 



(1) 



where, in our geometry (see Figure 
P± ( z ) = Pz {z) , 

P\\ (*) = (Px (Z)+ P y (z))/2. 

The dominant contributions to the integral in Eq. 
occur near the interface; in the bulk away from the in- 
terface, p± — p\\ and the integrand vanishes. For the 
specific case shown in Figure ^ where the interface sep- 
arates a bulk liquid from its corresponding vapor phase, 
the integral in Eq. can be replaced with an ensemble 
average of the difference between the normal and tangen- 
tial pressures, 



7p = -^(Pt-P\\) = y 



(P Z ) 



(Px) + (Py 



(2) 



The outer factor of 1 /2 in Eq. J5J) accounts for the pres- 
ence of two liquid- vapor interfaces. 

In most numerical simulations, the interatomic and 
electrostatic interactions are only applied within a cutoff 
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TABLE I: Parameters for commonly used three- and four-point models of water 



Parameter 


SPC/E 


TIP3P 


TIP3P-C 


TIP3P-Ew 


TIP4P 


TIP4P-Ew 


1H 


0.410 


0.417 


0.417 


0.415 


0.520 


0.5242 


qo 


-0.820 


-0.834 


-0.834 


-0.830 
















-1.040 


-1.0484 


Ohoh, deg 


109.47 


104.52 


104.52 


104.52 


104.52 


104.52 


Ioh, A 


1.0 


0.9572 


0.9572 


0.9572 


0.9572 


0.9572 












U. 1 JUU 


n 1 9^0 


eoo, kcal/mol 


0.1553 


0.1521 


0.1521 


0.102 


0.1550 


0.16275 


(TOO, A 


3.166 


3.1506 


3.1507 


3.188 


3.1536 


3.16435 


eoH, kcal/mol 






0.0836 








0"OH, A 






1.7753 








tHH, kcal/mol 






0.0460 








ohh, A 






0.4000 









range r c . The introduction of the cutoff in the interpar- 
ticle potential reduces the surface tension in much the 
same way that the introduction of a cutoff reduces the 
bulk pressure at constant density. Thus, the simulation 
result 7 P will underestimate the total surface tension; a 
better estimate of the total surface tension can be ob- 
tained from 



7 = 7p + Itail, 



(3) 



where TtoiK the tail correction for j p , can be determined 

from EuEa 



itail ' 2 



oo /•! />oo 



— oc J — 1 J r 



r 3 U' (r)g(r) (l - 3s 2 ) x 



p(z)p{z — sr) - (p G (z)) ) dr ds dz, (4) 



where U (r) is the pairwise potential, g (r) is the radial 
distribution function, p (z) is the observed interfacial pro- 
file, and pa (z) is a Gibbs dividing surface: 



PG (z) 



Ap 



sgn (z) . 



(5) 



Although the use of the tail correction in Eq. © im- 
proves the estimate of the surface tension, its use is re- 
stricted to systems in which the two phases contain the 
same components; for composite systems, such as water 
at the surface of a solid, only j p should be used. 

We assume in Eqs. (QJ and © that the liquid- vapor in- 
terface is centered at z = 0. In Eq. p c = (pi + p v ) /2 
is the average density of the two phases, and Ap = pi — p v 
is the difference between the average densities of the two 
phases. Thus, pc (z) — p v for z < and pi for z > 0. 
There are multiple possible choices for determining the 
observed density profile p(z). Although it is possible to 
use the profile calculated from the simulation directly, 
both the tail-correction and capillary-wave calculations 
are simplified by fitting the profile to a function. In the 
present study, the density profile is fitted to both an error 
function and a hyperbolic tangent function, as discussed 
in the following section. 



2. Capillary-wave method 

The thermodynamic approach for computing surface 
tension assumes a sharp liquid-vapor interface when in 
reality it is quite rough. The roughness of the inter- 
face increases at high temperatures, as seen in Figure ^ 
A second method for computing the surface tension as- 
sumes that the observed magnitude of the fluctuations is 
derived from two sources: an intrinsic contribution plus a 
logarithmic term that represents broadening of the inter- 
face as a result of the capillary waves 0,H,E3>0jEHE3- 

If the contributions from capillary waves can be decou- 
pled from density fluctuations, then the surface tension 
can be computed by determining the interfacial profiles 
for a number of different system sizes. The relationship 
between the observed interfacial width A and the intrin- 
sic interfacial width Aq is given by 



A 2 , 



k B T 



■In 



L 



Bo 



(6) 



where Lu = L x = L y is the length of the interface, and 
Bq is a characteristic length scale related to the short- 
wavelength cutoff in the interfacial behavior. It is un- 
necessary to determine Bq before computing the surface 
tension j w . 

Computation of j w requires the scaled density profile, 



9{z) 



PL - PV 



p(z) 



PL + PV 



(J) 



in the z-direction. Given ^(z), the variance in the 
derivative of the profile / (z) = iff' (z) can be computed 
in either real or reciprocal (Fourier) space: 



A 2 = 



JZoffW dz 



1 



/(0) 



d 2 f (q) 
dq 2 



(8) 



J 9=0 



where / (q) is the Fourier transform of f(z). The sim- 
ple form of Eq. © suggests that fitting the profile $ (z) 
to a functional form will be both more convenient and 
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lead to more accurate results than using the raw profile 
data. Several different functional forms for iff (z) have 
been proposed in the literature. Relying on mean-field 
arguments, most previous theoretical and computational 
studies of surface tension have fitted the profile to a hy- 
perbolic tangent function pllsillill . 



\ff t ( z ) = tanh 



(9) 



while Huang and Webb |43j and Beysens and Robert [44 
propose the use of an error function, 



* e (z) = erf 



(10) 



If the density profile \& (z) is fitted to a hyperbolic tan- 
ent function Eq. 0, then from Eq. JHJ we find that 



A 2 t = 7r 2 Wt 2 /48, 

while for an error function Eq. p0(l , the interfacial width 
Ag is given by 

A 2 e = w 2 e /27T. 

We will show that there is a significant discrepancy be- 
tween the surface tensions obtained from the hyperbolic 
tangent profile, Eq. J^Jl, and the error function profile, 
Eq. JTUJ), with the error function giving results in closer 
agreement with Eq. J2J). 



B. Water models 

We consider four different three-point models: the 
SPC/E model; the original TIP3P model; the modifica- 
tion of the TIP3P model H3 implemented in CHARMM 
(hereafter referred to as TIP3P-C); and the TIP3P-Ew 
model , a recent reparameterization incorporating the 
effects of Ewald summation. The parameters for the dif- 
ferent water models are summarized in Table [I] 

The basic structure of the different models are similar. 
The common features of all models include a specified 
oxygen-hydrogen bond length Ion and hydrogen-oxygen- 
hydrogen bond angle Ohoh, a charge on each hydrogen 
atom, and a Lennard- Jones 12-6 potential describing the 
interaction between the oxygen atoms, 



Ulj (roo) = 



^OO 











r- 




\roo 




\ roo J 




0, 





roo < r c 



roo > r c 
(11) 

where eoo and aoo are the model-dependent well depth 
and equilibrium 0-0 distance, and roo is the distance 
between oxygen atoms. The TIP3P-C model incorpo- 
rates hydrogen-hydrogen and hydrogen-oxygen Lennard- 
Jones interactions as well. 



In addition to the Lennard- Jones interaction, there are 
electrostatic interactions between the charge sites: 



U, 



N N 

i=i j=i 



(12) 



where q a is the charge on atom a, and is the distance 
between atoms i and j in the simulation box. Previ- 
ous studies have shown that significant variations in the 
values obtained for surface tension can occur depending 
upon how the sum in Eq. (|12f> is performed |3Sq. Except 
in Section Till CI Ewald summations were used through- 
out our simulations. 

We also consider a pair of four-point water models: the 
TIP4P model [13, and the TIP4P-Ew model 0, a re- 
cent reparameterization designed to account for the pres- 
ence of long-range interactions. The four-point models 
introduce a bare charge at a new site, designated M, lo- 
cated on the bisector of the HOH bond angle; the charge 
is of strength qM- The forces acting on the massless 
site are distributed to the O and H atoms in the same 
molecule |48| : 

Fy,o = (1 — 2a) Fij.M, 

where a = Iom / (Ioh cos (9hoh /2)) and Fy : M is the 
force acting on the massless site associated with oxygen i 
due to atom j 1 . For the TIP4P model, the charge is 
located Iom = 0.15 A away from the oxygen atom. The 
TIP4P-Ew model changes the values of Iom, the charge 
qM, as well as the separation aoo and well-depth eoo of 
the Lennard- Jones interaction. 



C. Simulation method 

1. Thermodynamic method 

To determine the surface tension of the various three- 
point water models, 1000 molecules were placed into a 
periodic, rectangular box of dimensions L x = L y = L\\ — 
2.3 nm and L z = L±_ = 13.5 nm. The increased system 
size in the z-direction minimizes the interactions of wa- 
ter molecules in the liquid phase with their z-periodic 
images through the long-range Coulombic interactions 
in Eq. (|12fl . Similarly, 1000 molecules of the four- 
point models were simulated in a box with dimensions 
L x = L y = L|| = 2.7 nm and L z = Lj_ = 12.0 nm, 
each also containing 1000 molecules. The initial config- 
uration was constructed by placing the water molecules 



1 In the case of four-point water models, while the contribution 
to the virial from the Lennard-Jones interactions can be com- 
puted using the more computationally efficient form X]i r i ' 
the contribution to the virial from short-range electrostatic forces 
is computed using th.6 palrwisc form jy-j 



at the center of a simple cubic lattice with 7 molecules 
each in the x- and y-directions, and the z-spacing cho- 
sen to create a density of 0.98 gem -3 for the three-point 
models, and 1.00 gem -3 for the four-point models. The 
same starting configuration was used for all simulations 
of a given water model. At equilibrium, the thickness of 
the slab in the z-direction varied between approximately 
5.5 nm at 300 K and 7.5 nm at 500 K. 

For each of the six models examined, molecular dynam- 
ics (MD) simulations were performed in the NVT en- 
semble in 25-degree increments between 300 K and 500 K 
using the LAMMPS simulation package |50|. The cut- 
off for the Lennard- Jones potentials and the short-range 
cutoff for the electrostatic potentials were set to 10 A, 
unless otherwise specified. The bond lengths and bond 
angles of the various models were constrained using the 
SHAKE technique [HlJ . The equations of motion were in- 
tegrated using the Verlet algorithm with velocity rescal- 
ing to control the temperature. The difference in the 
surface tension between simulations performed with ve- 
locity rescaling and those with a Nose-Hoover thermostat 
with a damping constant of 100 ps _1 was significantly less 
than the simulation uncertainty. Each simulation was 
performed for a total of 2 ns with time step At = 1 fs . 
The system was allowed to equilibrate for 1 ns; data from 
the second 1 ns were used to compute the surface tension. 

The electrostatic interactions were calculated using 
the particle-particle particle-mesh (PPPM) technique of 
Hockney and Eastwood |5§j. The mesh spacing in this 
work was selected to ensure that the root-mean-squared 
accuracy of the force calculation was within 10~ 4 ; the 
resulting grid was of dimensions 12 x 12 x 48. Most 
previous simulations were carried out with a maximum 
of hf ax = 20 cells in the z-direction HJ H3- In sev- 
eral of those studies, simulations were carried out with 
/i™ = 10 or less, and some did not include long-range 
electrostatic interactions at all [3j| . We consider the ef- 
fects of mesh refinement on the surface tension in Section 




-2 2 

z position (nm) 

FIG. 2: Density profiles for the SPC/E model of water at 300 
K (thick solid line), 400 K (thick dashed line), and 500 K 
(thick dashed-dotted line). Fits to error functions are shown 
as thin solid lines. 



last 750 ps of the simulation were used for recording data; 
the preceding steps were used for equilibration and dis- 
carded. We output the position of every atom after every 
5000 timesteps, and then assigned each atom to one of 
500 bins depending on its location in the z-direction. To 
ensure that the interfacial profile was not altered by drifts 
in the density profile, the profile was shifted so that the 
center of the mass was located at z = 0. Sample profiles 
for the SPC/E model of water at several temperatures 
are shown in Figure El After the average density profile 
p (z) was computed, the two halves of the profile, on ei- 
ther side of z = 0, were averaged together, rescaled to 
values between —1 and 1 using Eq. JJJ), and then fitted 
to hyperbolic tangent and error functions of the form of 
Eq. © and Eq. JTOJ. 



III. THERMODYNAMIC SURFACE TENSION: 
RESULTS AND DISCUSSION 



2. Capillary-wave method 



A. Temperature dependence 



Observation of capillary waves requires simulations 
with larger interfacial surface areas than were used in 
the thermodynamic method above. Consequently, we 
studied systems with L x — L y — Lu varying between 
9.2 nm and 46.0 nm. The resulting simulations used to 
compute the surface tension have surface areas between 
84.6 nm 2 and 2116.0 nm 2 , and contained between 16,000 
and 400, 000 water molecules. To construct the initial 
configuration, we take an equilibrated sample and repli- 
cate it multiple times in the x- and y-directions. The 
SPC/E model was used for this study, as it was the most 
computationally efficient of the models studied. 

To ensure that artifacts from the replication process 
were eliminated, the simulation time varied between 
1.0 ns and 6.0 ns as an increasing function of Lu . Only the 



The surface tension 7 = 7 p + jtaii of the various water 
models, computed using Eq. (J2J), are collected in Table 
ITT1 Results for the three-point models are shown in Fig- 
ure^ a) as a function of temperature. Using the method 
of Flyvbjerg and Petersen [53, the uncertainty in the re- 
sults was found to be between 2.4 and 3.0 mN / m. Com- 
paring the four three-point models, we find that SPC/E 
model is the closest to the experimental data, with bet- 
ter agreement at higher temperatures. For most temper- 
atures, the various TIP3P models agree with the SPC/E 
model within the uncertainty of the simulation. Al- 
though the three-point models considered do not achieve 
agreement with experimental data, the overall tempera- 
ture dependence of 7 for the models is in good agreement 
with the experimental data. For the three-point mod- 
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TABLE II: Surface tension" for three- and four-point water models, including tail correction 



Surface Tension, 7 (niN/m) 



TpTTl T"_PT"£li~ 11 VP ( T\ 1 
_I_Ci.ll 0, L Lll " \ IV 1 


SPC/E 


TIP3P 


1 11 Ol _L_J W 


TIP3P-C 


TIP4P 


TTP4P-F,w 

_L 11 4:1 _____ W 


Expt. 


oUU 


00.4 


01.1 


A7 A 


A ft 8 


Oo.D 


Dl.Z 


71 7 


325 


47.9 


45.9 


43.6 


47.5 


50.9 


55.6 


67.6 


350 


47.0 


42.8 


39.2 


45.7 


45.7 


52.7 


63.2 


375 


44.6 


37.8 


37.7 


39.9 


41.4 


48.1 


58.4 


400 


37.6 


35.5 


34.3 


36.9 


35.9 


43.5 


53.3 


425 


32.0 


31.5 


28.9 


31.9 


31.2 


38.6 


47.9 


450 


30.6 


27.3 


25.9 


28.2 


25.7 


34.5 


42.1 


475 


26.8 


24.7 


23.2 


23.2 


19.1 


29.3 


36.0 


500 


23.2 


17.0 


16.9 


18.1 


15.2 


24.8 


29.5 



Notes: "Uncertainty for all simulation results is between 2.4 and 3.0 mN/m. b Experimental data taken from Refs. 
J2a.l27-,l2ff : data above 400 K is extrapolated from quadratic fit provided in Ref. |27|I . 



TABLE III: Interfacial properties of SPC/E water as a function of temperature, for r c — 10 Aand Ly = 2.3 nm 



Temperature 
T, K 


Interface thickness 

w e , A 


Liquid density 
PL, g/cm 3 


Vapor density 
Pv, g/cm 3 


Tail correction 
jtaii, mN/m 


300 


3.12 


0.990 


0.0005 


5.5 


325 


3.37 


0.977 


0.0008 


5.2 


350 


3.75 


0.959 


0.0005 


5.0 


375 


4.22 


0.941 


0.0005 


4.8 


400 


4.63 


0.913 


0.0015 


4.5 


425 


5.23 


0.886 


0.0023 


4.1 


450 


5.74 


0.855 


0.0049 


3.8 


475 


6.00 


0.818 


0.0093 


3.5 


500 


7.54 


0.779 


0.0199 


3.0 



els, 7 is generally between 15 mN / m and 20 mN / m less 
than the experimental data, especially in the tempera- 
ture range 300 K < T < 425 K. 

Alejandre et al. j^H and Shi et al. reported excel- 
lent agreement with experimental results for the SPC/E 
model. However, our results for the SPC/E model clearly 
disagree with their data as well as with experimental re- 
sults, although the simulations were performed under es- 
sentially identical conditions with respect to the number 
of molecules and the dimensions of the system, potentials 
employed, temperature range, and cutoffs. We study the 
potential causes of the disagreement in the results below 
in Sections ImCl and InTDl 

Results for the four-point models as a function of tem- 



TABLE IV: Surface tensions j p and 7 and liquid-phase den- 
sity pL for the SPC/E model as a function of LJ cutoff r c , 
with and without tail correction, at 300 K 



r c (A) 


7 P (mN/m) 


7 (mN/m) 


p L (g/cm 3 ) 


10.0 


46.3 


51.8 


0.990 


12.0 


51.2 


55.0 


0.992 


14.0 


47.9 


50.6 


0.994 


16.0 


49.7 


51.8 


0.996 


18.0 


49.9 


51.5 


0.996 


20.0 


52.8 


54.1 


0.995 



perature are shown in Figure IHIb). The uncertainty for 
the four-point models is the same as for the three-point 
models. Like the three-point models, the four-point mod- 
els underestimate the surface tension, with the TIP4P 
model offering results comparable to the TIP3P-C and 
TIP3P-Ew models, while the performance of the TIP4P- 
Ew model is significantly closer to the experimental data 
than any of the other models examined here. Unlike the 
TIP3P models, at most temperatures considered here, 
the TIP4P and TIP4P-Ew models do not agree within 
simulation uncertainty. 



TABLE V: Surface tension 7 P and 7 for the SPC/E model as 
a function of LJ cutoff 



7 P (mN/m) 7 (mN/m) 



T(K) 


r c = 10 A 


r c = 16 A 


r c = 10 A 


r c = 16 A 


300 


49.9 


50.2 


55.4 


52.3 


325 


42.8 


48.7 


47.8 


50.7 


350 


42.2 


47.2 


47.0 


49.1 


375 


40.0 


40.3 


44.6 


42.1 


400 


33.3 


37.7 


37.6 


39.4 


425 


28.1 


33.5 


32.0 


35.0 


450 


27.0 


29.2 


30.6 


30.6 


475 


23.6 


23.6 


26.8 


24.9 


500 


20.5 


22.3 


23.2 


23.3 
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can evaluate Eq. After integration over s, 
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FIG. 3: (Top figure) Surface tension of the three-point models 
of water as a function of temperature: SPC/E (filled circles), 
TIP3P (squares), TIP3P-C (triangles), and TIP3P-Ew (dia- 
monds). Simulation data from Alejandre el al. [3111 (open 
circles) are included for comparison. (Bottom figure) Surface 
tension of the four-point models of water as a function of tem- 
perature: TIP4P (circles) and TIP4P-Ew (squares). In both 
figures, experimental data jg^, [2^ (solid curve) and ex- 
trapolation of quadratic fit to higher temperatures (dashed 
curve) are included for comparison. 



B. Tail correction 



While tail corrections can exist for both the Lennard- 
Jones and the electrostatic interactions, by using Ewald 
summations, we avoid the need for a Coulombic tail 
correction. In evaluating the tail correction, Eq. Q 
for the Lennard- Jones potential in the region r a, 
dU (r) jdr w 24ecr 6 r~ 7 . Assuming that for r > r c , the 
radial distribution function g (r) w 1 and that the den- 
sity profile is an error function of the form Eq. (|10|l , we 



Jtail 

12nea e 



oo /»oo 



erf (gz) p -e{r+z)\ 



3 r 7 



(^z (r-z) + e 4 « 2 " (£ 2 z (z + r) + l) - l) 

- Z -^T 1 ( 3 -e{z 2 -r 2 ))f (f, r, z) dr dz, (13) 

where £ = y/ir/w e and 

/ (£, r, z) = erf (£ (r + z)) + erf (£ (r - z)) . 

Eq. \i'6\ is then evaluated numerically for each den- 
sity profile. For the SPC/E model of water, the interfa- 
cial thickness, density difference, and tail correction are 
shown in Table ITTT1 Because of the existence of capillary 
waves at the interface, as shown in Section IIVI the in- 
terfacial thickness parameter w e depends logarithmically 
on the length Lu of the interfacial surface. 

The tail correction depends strongly upon the 
Lennard- Jones parameters, scaling as eer 6 , and decays 
exponentially as a function of the chosen interaction cut- 
off r c . Because eoo an d a oo are approximately equal 
for the SPC/E and TIP3P models, the tail corrections 
at all temperatures are almost identical for the two mod- 
els. While the TIP3P-C has additional tail corrections 
for the O-H and H-H interactions, their magnitudes are 
negligible in comparison to the correction for the 0-0 
interaction. Only the TIP3P-Ew model, which has a 
significantly smaller value for the Lennard-Jones inter- 
action strength £oo, has a noticeably different tail cor- 
rection from the other three-point models. The tail cor- 
rections for the four-point models are likewise close to 
that of the SPC/E model, with the TIP4P model hav- 
ing a slightly smaller tail correction and the TIP4P-Ew 
model a slightly larger tail correction. 



C. Cutoff effects 

The original parameterizations for the SPC/E and 
TIP3P models of water employed cutoffs for both electro- 
static and Lcnnard-Jones interactions Elini. To study 
the effect of varying the electrostatic cutoff, we applied 
a short-range cutoff to both the LJ and electrostatic po- 
tentials of the SPC/E model, truncating the potentials at 
r cut = 10, 12, 14, 16, 18, and 20 A. Using Eq. © to de- 
termine the surface tension, we found that the estimated 
values of the surface tension were nonsensical, ranging 
between —3700 and 1700 mNm -1 , with no value smaller 
in magnitude than 140 mNm -1 . This shows that trun- 
cated electrostatic potentials are inappropriate for use in 
the determination of the surface tension of water. 

To determine the effect of varying only the range of the 
Lennard-Jones interaction on 7 P before incorporating the 
tail correction, we performed runs for the SPC/E model 
at 300 K with LJ cutoffs of 10, 12, 14, 16, 18, and 20 A, 



using the PPPM Ewald technique for the electrostatic 
forces. The starting configuration for these runs was 
the final configuration from the 300 K simulation used 
to compute the surface tension in Section IIII Al The 
resulting values of the surface tension are shown in Ta- 
ble IIVI The values for the surface tension are within 
the simulation uncertainties, although they tend to rise 
with increasing r c . This is further reflected in the den- 
sity profiles, which show that the liquid-phase density Pl 
increases with r c for values of r c between 10 A and 18 A. 
The overall effect of the choice of r c can be seen when 
comparing the surface tensions of SPC/E water in the 
temperature range 300 K to 500 K for LJ cutoffs of 10 A 
and 16 A. The resulting data are shown in Table W\ The 
data demonstrate that the surface tensions for 10 A and 
16 A cutoffs are equal, within simulation uncertainty, af- 
ter the corresponding tail corrections have been applied 
to each set of data. Thus, the 10 A Lennard- Jones cutoff 
with long-range tail corrections is sufficiently accurate for 
computing the surface tension. 



D. Reciprocal space accuracy dependence 

Alejandre et al. |3l| found that the surface tension de- 
pended on the mesh refinement |h max | used to evaluate 
the long-range Coulombic interactions. To test the de- 
pendence of 7 P on |h max | , we show results in Figure 01 
for three different models as a function of |h max |. From 
the figure, several trends become apparent. First, for 
both the TIP3P and TIP3P-Ew models, the long-time 
average of the surface tension depends significantly on 
the value |h max |: the long-time average for a 5 x 5 x 20 
mesh (|h max | = 20, rms accuracy 4.0 x 10~ 3 ) is between 
15 and 20mNm~ 1 larger than for a 12 x 12 x 48 mesh 
(|h max | = 50, rms accuracy 10 -4 ). For larger values of 
|h max | , there is no significant adjustment in the surface 
tension. Additionally, as found above, for a given value 
of the precision, there is little difference in the long-time 
average of the two models. 

For the SPC/E model, there is little difference between 
the equilibrium values for the 5 x 5 x 20 mesh and the 
12 x 12 x 48 mesh. However, it takes approximately 1 ns to 
achieve agreement between the two precision levels; be- 
fore this, the less-refined mesh has a significantly greater 
surface tension. In the work of Alejandre et al., the total 
simulation time was only 0.375 ns . After 0.375 ns, the 
average surface tension from our simulation using the 
5x5x20 mesh was approximately 7 = 60 mN m~ 1 , which 
closely corresponds with the result 7 = 60.6 mNm -1 ob- 
tained in the earlier study. However, for long simulation 
times, the surface tension of the SPC/E model does not 
exhibit a strong dependence on the mesh size. Conse- 
quently, since the long-time averages are essentially equal 
within simulation uncertainty, we have used the finer 
12 x 12 x 48 mesh refinement for all of the simulations 
reported in this paper unless otherwise specified. 

To further demonstrate that the accuracy of the 
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FIG. 4: (Color online) Equilibration of the surface tension of 
the SPC/E (black circles), TIP3P (red squares), and TIP3P- 
Ew (blue diamonds) models of water for rms k-space accu- 
racies of 4 x 1(T 3 (dashed curves, open symbols) and 10 
(solid curves, solid symbols). 




Temperature (K) 



FIG. 5: Surface tension for the SPC/E model of water as a 
function of mesh accuracy averaged over (top) 125-375 ps and 
(bottom) 1-2 ns for a 5 x 5 x 20 mesh (circles) , a 12 x 12 x 48 
mesh (diamonds), and a 20 x 20 x 80 mesh (triangles). 



reciprocal-space calculation does not affect the results 
for the surface tension of the SPC/E model when aver- 
aged over sufficiently long times, we performed simula- 
tions of the SPC/E model, as described above, with k- 
space meshes of 5 x 5 x 20, 12 x 12 x 48, and 20 x 20 x 80, 
corresponding to rms accuracies of 4 x 10~ 3 , 10~ 4 , and 
10 -5 , respectively. The results are shown in Figure [S] for 
the time ranges of 125 to 375 ps, used by Alejandre et 
al, and 1 to 2 ns, used in the present work. Clearly, the 



TABLE VI: Comparison of surface tension for the SPC/E 
model as a function of calculation method 

Surface tension, mN/m 
Capillary- wave 

T (K) Erf, y w Tanh, j w Pressure-integration, 7 P 



300 46.1 36.5 49.9 

400 32.0 27.8 33.1 

500 19.0 19.0 20.2 



TABLE VII: Surface tension for the SPC/E model as a func- 
tion of system size at 300 K 

Interfacial length Ly, A j p (mN/m) 7 (mN/m) 

2.3 49l) 55.4 

9.2 52.0 57.5 

11.5 51.5 57.0 

13.8 51.8 57.3 

16.1 51.9 57.4 

23.0 51.6 57.1 

34.5 52.0 57.5 

46.0 51.8 57.3 



results for the shorter time range (125 to 375 ps) have not 
converged: there is significant disagreement of as much 
as 20 mN / m, particularly at higher temperatures. How- 
ever, for longer times, the results have converged, with 
the differences among the three mesh refinements being 
essentially within simulation uncertainty. It is interest- 
ing to note that the converged surface tensions of the 
different mesh refinements corresponds very closely to 
the profile obtained for the intermediate mesh refinement 
(12 x 12 x 48) after 375 ps. Our intermediate mesh re- 
finement is comparable to the most refined k-space mesh 
considered by Alejandre et al. [3lJ (shown in Fig. 7 of 
their paper for a single, instantaneous configuration). 



IV. CAPILLARY WAVES 

The interfacial width A 2 was computed for various sys- 
tem sizes as described in Section [H] The resulting plots 
of A 2 versus In L were computed and the value of *y w ex- 
tracted using Eq. JJjJ. A plot showing the data obtained 
for 300 K for the SPC/E water model is shown in Figure 
|HJ the resulting values of 7^ for the two functional forms 
at 300 K, 400 K, and 500 K are shown in Table IVD 

Our results indicate that at all three temperatures, 
good agreement between the pressure-integration esti- 
mate of the surface tension, j p , and the capillary- wave 
estimate of the surface tension, j w , is obtained only if 
the interfacial density profile is fit to an error function. 
There is substantial disagreement between the error func- 
tion and hyperbolic tangent functions at lower tempera- 
tures: the hyperbolic tangent profile always yields lower 
estimates than the error function profile. However, the 
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FIG. 6: Regression fit of Eq. JHj for hyperbolic tangent (dia- 
monds) and error function (circles) profiles at 300 K. 



magnitude of the discrepancy between the two estimates 
of 7 decreases the two narrows as temperature increases, 
and essentially vanishes at 500 K. 

In addition to computing the surface tension j w , we 
can also compute an upper-bound estimate for the in- 
trinsic interfacial width Ao |8j . After obtaining the slope 
a w and intercept f3 w for the least-squares fit of Eq. J(|J|, 
we assume that the parameter Bq — cAo, where c is a 
constant to be specified. The intercept (3 W and A are 
then related by 

f3 w = A 2 - a w In (cA ) . (14) 

For values of c less than some threshold c* , there is no real 
solution to ^14(1: above the threshold, Ao quickly decays 
as c increases. Thus, Ao has a maximum at the threshold 
value c = c* where the imaginary part of the solution 
vanishes. For the system sizes under consideration, we 
find that the maximum intrinsic width Ao ~ 0.8 A, 1.0 A, 
and 1.5 A at T = 300 K, 400 K, and 500 K. 

Further evidence that fitting the density profile to an 
error function yields more accurate results than fitting 
to a hyperbolic tangent function can be seen by compar- 
ing the density fits themselves. As shown in Figure [7| 
although the two profiles are similar, the error function 
fit more closely adheres to simulation results than the 
hyperbolic tangent profile. Although the x 2 parameter 
for both functional fits was relatively small, the coeffi- 
cient for the hyperbolic tangent x 2 parameter (~ 10~ 3 ) 
was approximately two orders of magnitude larger than 
the corresponding error function x 2 parameter (~ 10 -5 ). 
The greater accuracy of the error function is further seen 
by comparing the magnitude of the differences between 
the simulation results and the fitted functional profiles, 
as shown in Figure |H1 

In addition to observing capillary wave behavior, the 
larger simulations can be used to study the effect of the 
interfacial area L 2 on the surface tension computed using 
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FIG. 7: Density profile in the interfacial region for a slab 
of 400,000 SPC/E water molecules at 300 K (L\\ = 46.0 nm). 
Simulation results are shown as circles; fits to error function 
(solid curve) and hyperbolic tangent function (dashed curve) 
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we expect the simulation values for the surface tension 
to remain constant, independent of L». Examining the 
results for the surface tension 7 P (without tail correction) 
versus the system size N, as shown in Table IVTT1 we find 
a slight increase in 7 as Ly increases, although the results 
remain within simulation uncertainty, even for large val- 
ues of Lu. 

V. CONCLUSIONS 

An investigation of surface tension as a function of tem- 
perature for a number of popular three- and four-point 
water models shows systematic disagreement between ex- 
perimental and simulation results for all six models con- 
sidered. All six models considered consistently underes- 
timate the surface tension relative to the experimental 
data. The TIP4P-Ew model is closest to the experimen- 
tal data, although not in quantitative agreement. 

Examining the SPC/E model in greater detail, we have 
illustrated the importance of having a sufficiently long 
simulation time and sufficiently fine k-space mesh: sig- 
nificant variations in the surface tension can result if the 
equilibration period is too short. In the SPC/E, TIP3P, 
and TIP3P-Ew models, significant variations were also 
observed if a small number of k-space vectors are used. 

The present study also includes the first in-depth study 
of the effect of large system sizes on interfacial proper- 
ties, studying capillary waves for systems of up to 4 x 10 5 
molecules. We have demonstrated that examining cap- 
illary waves at the liquid-vapor interface can be used to 
determine the surface tension of real fluids, and that the 
use of an error function profile offers better results in 
comparison to experimental data than the use of a hy- 
perbolic tangent profile. Finally, we note that the inter- 
facial width A depends on the interfacial length Lm, and 
that attempting to extract an interfacial width without 
taking into account the effect of capillary waves [2j| is 
incorrect. 



FIG. 8: Difference between simulation results and the er- 
ror function (solid curve) and the hyperbolic tangent (dashed 
curve) fits for the same system shown in Figure 
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